

//@#ifndef num_threads
//    @#define num_threads = 1
//@#endif

// Disable parallel processing explicitly
//options_.parallel_info.number_of_threads = 1;
//options_.parallel = false;


var w1,w2,T1,T2,pi12,pi21,x1,x2,r1,r2,taux,pi11,pi22,P2,riskfree1,riskfree2,rmul;

varexo sa1,sa2;

parameters betad, theta1, dbar, sigmaH, alpha2, L1, L2,delta1, rhoa, rho,alpha1,gL,std_ea,rbest,extax_ss,lrss,Ts1,Ts2,xss1,xss2
            Pss2,wss1,wss2,rss,piss12,piss21,riskfree_ss,sigmaF;

load setpara_ramsey.mat

ncountries     =  mpara(1); % number of countries
nsectors       =  mpara(2);
betad   = mpara(3);     % discount factor
theta1  = mpara(4);
sigmaH = mpara(5);
gL      = mpara(6);
dbar   = mpara(7);
delta1 = mpara(8);
L1     = mpara(9);
L2      = mpara(10);
alpha1 = mpara(11);
alpha2  = mpara(12);
rhoa    = mpara(13);
std_ea  = mpara(14);
extax_ss    = mpara(15);
lrss   = mpara(16);
Ts1    = mpara(17);
Ts2   = mpara(18);

xss1 = mpara(19);
xss2 = mpara(20);
Pss2 = mpara(21);
wss1 = mpara(22);
wss2 = mpara(23);
rss = mpara(24);
piss12 = mpara(25);
piss21 = mpara(26);
riskfree_ss = mpara(27);
sigmaF = mpara(29);



rho   = betad*(1-delta1);
rbest = delta1/(theta1*(1-betad*(1-delta1))+delta1);



model;

[name = 'P1 definition']
(exp(T1)*exp(w1)^(-theta1)+exp(T2)*(exp(w2)*dbar)^(-theta1))^(-1/theta1)=1;

[name = 'P2 definition']
P2 = -1/theta1*log(exp(T1)*(exp(w1)*(1+taux)*dbar)^(-theta1)+exp(T2)*(exp(w2))^(-theta1));

[name = 'balance of payment']
pi12*exp(x1)=pi21*exp(x2);

[name = 'optimal innovation']
exp(w2)/(exp(sa2)*alpha2) = 1/theta1*exp(w2)*(L2-exp(r2)*L2)/exp(T2)+betad*(1-delta1)*(exp(x2(+1))^(-sigmaF)*exp(P2(+1))^(sigmaF-1)*exp(w2(+1))/(exp(sa2(+1))*alpha2)*exp(x2)^sigmaF*exp(P2)^(1-sigmaF));

[name = 'expenditure x2']
exp(x2) = (1+theta1)/theta1*exp(w2)*(1-exp(r2))*L2;

[name = 'expenditure x1']
exp(x1) = (1+theta1)/theta1*exp(w1)*(1-exp(r1))*L1+taux/(1+taux)*pi21*exp(x2);

[name = 'evolution T1']
exp(T1) = exp(T1(-1))*(1-delta1)+exp(sa1)*alpha1*exp(r1)*L1;

[name = 'expenditure T2']
exp(T2) = exp(T2(-1))*(1-delta1)+exp(sa2)*alpha2*exp(r2)*L2;


//pi11 = exp(T1(-1))*(exp(w1))^(-theta1)/(exp(T1(-1))*exp(w1)^(-theta1)+exp(T2(-1))*(exp(w2)*dbar)^(-theta1));
//pi22 = exp(T2(-1))*(exp(w2))^(-theta1)/(exp(T1(-1))*(exp(w1)*(1+taux)*dbar)^(-theta1)+exp(T2(-1))*(exp(w2))^(-theta1));

pi11 = exp(T1)*(exp(w1))^(-theta1)/(exp(T1)*exp(w1)^(-theta1)+exp(T2)*(exp(w2)*dbar)^(-theta1));
pi22 = exp(T2)*(exp(w2))^(-theta1)/(exp(T1)*(exp(w1)*(1+taux)*dbar)^(-theta1)+exp(T2)*(exp(w2))^(-theta1));
pi12 = 1-pi11;
pi21 = 1-pi22;


riskfree1  = exp(x1)^(-sigmaH)/(betad*exp(x1(+1))^(-sigmaH));
riskfree2  = exp(x2)^(-sigmaF)*exp(P2)^(1-sigmaF)/(betad*exp(x2(+1))^(-sigmaF)*exp(P2(+1))^(1-sigmaF));
  
//sa1 = rhoa*sa1(-1)+std_ea*ea1;
//sa2 = rhoa*sa2(-1)+std_ea*ea2;

rmul = (1+theta1*pi22)/(1+taux)-theta1*pi22;

end;


planner_objective(exp(x1)^(1-sigmaH)/(1-sigmaH));

ramsey_model(instruments=(taux,r1),planner_discount=betad); 


initval;

T1 = Ts1;
T2 = Ts2;

w1 = log(wss1);
w2 = log(wss2);
pi12 = piss12;
pi21 = piss21;
x1 = log(xss1);
x2 = log(xss2);
r1 = log(rss);
r2 = log(rss);
taux = extax_ss;
pi11 = 1-piss12;
pi22= 1-piss21;
P2 = log(Pss2);
riskfree1 = log(riskfree_ss);
riskfree2 = log(riskfree_ss);
rmul = 0;

sa1 = 0;
sa2 = 0;
end;

endval;

sa1 = 0;
sa2 = 0;

taux = extax_ss;
r1 = lrss;

end;
steady;
//resid; 


perfect_foresight_setup(periods=2000);
perfect_foresight_solver(stack_solve_algo=0);

/*
rplot x1;
rplot T1;
rplot T2;
rplot taux;
rplot riskfree1;
rplot riskfree2;
rplot x2;
rplot sa1;
rplot w1;
rplot w2;
rplot P2;
*/

//stoch_simul(order=2, periods = 1000, irf = 0);
//stoch_simul(order=2,irf=40) sa1,T1,T2,taux,r1,r2,x1,x2,riskfree1,riskfree2,w1,w2,P2;
//evaluate_planner_objective;

//check;

//write_latex_dynamic_model;
//write_latex_static_model;
//write_latex_definitions;
//write_latex_parameter_table;

//collect_latex_files;

//model_diagnostics(M_,options_,oo_);


//rplot x1,T1,T2,taux,riskfree1,riskfree2,x2,sa1,pi11,pi22;



